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ABSTRACT 


iii 


A  computational  investigation  of  the  fracture  mechanics  of  metals  and  metal 
matrix  composites  has  been  carried  out.  The  ductile  fracture  of  structural  alloys  was 
analyzed  through  a  two  dimensional  non-linear  finite  element  approach,  while  the 
mechanics  of  load  transfer  in  silicon  carbide  (SCS-6)  fiber  reinforced  titanium  alloy  (Ti- 
15V-3Cr-3Al-3Sn)  were  studied  using  a  local-global  finite  element  analysis  procedure. 

The  computed  values  of  the  7-Integral  for  compact  tension  specimens  of  steel 
and  aluminum  alloys  (0.533  ^  dV/  ^  0.884)  remain  path  independent  up  to  a  certain 
load  which  is  attributed  to  crack  initiation,  and  then  diverge.  There  is  a  unique  signature 
of  the  strain  energy  density  (dWIdV)  ahead  of  the  crack  in  the  0*  direction:  the  strain 
energy  first  decreases,  reaches  a  minimum  and  then  increases  with  increasing  distance 
from  the  crack  tip.  The  minimum  strain  energy  shows  a  unique  dependence  on  the 
apfdied  load.  This  learLi  the  prediction  of  the  fracture  loads  for  the  cracked 
specimens. 

A  novel  coordinate  system  rotation  w'as  employed  in  extracting  the  boundary 
conditions  from  the  two-dimensional  global  model  to  the  three-dimensicHial  local  model 
for  the  local-global  finite  element  analv^is  of  the  unidirectional  composite.  The 
mechanics  of  load  transfer  and  the  subsequent  fracture  and  were  studied  over  a  range  of 
interface  strengths  ( 138  -  934  MPa).  The  analysis  leads  to  the  prediction  of  debonding 
betu'een  fiber  and  matrix.  The  degree  of  load  uansfer  from  the  matrix  to  the  fiber  was 
found  to  depend  upon  both  the  angular  orientation  and  distance  frmn  the  crack  front  An 
optimum  value  of  the  interface  strength,  (cfi)opi.  has  been  established  for  the  composite 
for  its  maximum  load  bearing  capability. 
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CHAPTER  1 


INTRODUCTION  TO  MICROMECHANISMS 

OF  FRACTURE  IN  COMPOSITES 


The  micromechanistns  of  fracture  in  polymeric,  metallic  or  ceramic  composite 
materials  reinforced  with  strong,  stiff  fibers  interact  in  a  complex  manner.  These 
include,  but  are  not  limited  to,  matrix  and  fiber  cracking,  crack  bridging,  delamination 
and  crack  deflection.  Metal  matrix  composites  present  the  additional  problem  of  yielding 
and  plastic  deformation  in  localized  regions.  Depending  on  the  composite  system  under 
investigation,  the  interface  will  affect  the  stress  distribution  around  a  loaded  crack  tip. 

In  the  past,  work  has  centered  on  experimental  observations  of  the  micromechanical 
causes  of  damage  [1*2],  closed  form  analytical  procedures  designed  to  calculate  laminate 
mechanical  properties  and  damage  accumulation  [3-8],  characterization  of  the  interphase 
region  through  finite  element  analysis  and  complex  constitutive  relations  [9-14],  and 
local-global  finite  element  calculations  of  cracked  composites  [15-17]. 

A  broad  base  of  knowledge  is  required  to  analyze  damage  in  composite 
materials.  Some  background  on  elastic-plastic  fracture  mechanics  will  first  be  given  to 
comprehend  how  plasticity  plays  a  role  on  both  the  microscopic  and  macroscopic  levels. 
Next,  the  metallurgical  and  composite  damage  mechanisms  in  these  composites  will  be 
categorized  and  quantified  to  account  for  the  multiple  energy  absorbing  mechanisms. 
Finally,  a  state  of  the  art  review  w  ill  be  given  on  analytical  and  computational 
micromechanical  modeling  techniques. 


Elastic-Plastic  Analysis  of  Fracture 


As  noted  by  Orowan  [18],  the  work  required  to  fracture  ductile,  tough  metal  is 
greatly  increased  due  to  an  additional  plastic  deformation  work  (yp).  This  deformation  is 
absorbed  through  a  process  of  microvoid  nucleation,  growth  and  linking. 
Microscopically  these  voids  may  begin  to  form  in  individual  grains  at  stresses  much 
lower  than  the  yield  point  E)epending  the  strain  hardening  s\nd  strain  rate  hardening  of  a 
material,  void  nucleation  and  growth  may  occur  homogeneously  in  areas  of  high 
dislocation  densities  or  may  occur  heterogeneously  along  grain  boundaries  or  particle 
inclusions  [19].  Conversely,  the  plastic  deformation  of  the  metal  on  a  large  scale  can  be 
experimentally  observed  on  pre-polished  and  loaded  samples  as  variations  in  the 
material's  surface  roughness.  Lee  [20]  found  that  the  intensity  of  reflected  light  is 
linearly  related  to  the  surface  roughness,  which  is  directly  proportional  to  the  effective 
strain.  This  allows  for  an  accurate  measurement  of  the  crack  tip  plastic  zone  size. 
Numerical  analyses  have  been  performed  verifying  the  extent  of  plastic  deformation  on 
globally  isotropic  material  models.  Criteria  such  as  the  7-Integral  are  then  used  to 
characterize  the  crack  tip  process  zone. 

Rice's  7-Integral  [21]  has  been  formulated  as  a  possible  fracture  criterion  for 
materials  that  exhibit  nonlinear  stress-strain  behavior.  7  has  been  shown  to  accurately 
predict  crack  initiation  in  structural  alloys  [22,23]  and  has  been  studied  numerically  as  it 
^plies  to  limited  crack  extension  {dJIda)  [24-28].  The  need  for  a  ductile  fracture 
criterion  stems  from  the  irtherent  plasticity  that  surrounds  a  crack  tip  prior  and  during 
crack  growth,  thus  exceeding  the  approximations  used  in  linear  elastic  fracture 
mechanics  (LEFM).  Elastic-plastic  fracture  mechanics  (EPFM)  must  be  employed  when 
the  crack  tip  plastic  zone  (CTPZ)  is  large  compared  to  the  dimensions  of  the  specimen. 
Fracture  parameters  such  as  the  crack  tip  op^ening  displacement  measurement  (CTOD) 
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and  analytical  solutions  of  the  slip-line  field  have  been  studied  as  they  relate  to  ductile 
instability  and  have  met  with  some  success  [29].  From  the  work  of  Hutchinson  [30], 
Rice  and  Rosengren  [31],  the  stresses,  strains,  and  displacements  are  known  within  a 
radius,  from  the  crack  tip  in  the  case  of  large  scale  and  general  yielding  for  linear 
hardening  elastic-plastic  materials.  This  zone  is  the  so-called  HRR  zone  (taken  from  the 
work  of  Hutchinson,  Rice  and  Rosengren),  while  the  material  outside  the  crack  tip 
plastic  zone  remains  in  the  elastic  region.  The  stresses  and  strains  within  the  HRR  zone 
may  be  described  by  the  following  equations; 
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The  constants  aocome  from  the  Ramberg-Osgood  description  of  a  material's 
stress  strain  curve  described  by  the  following  relationship: 


(1.4) 


The  work  of  Hutchinson,  Rice  and  Rosengren  gained  recognition  due  to  the 
variety  of  problems  that  could  be  described  by  the  so-called  HRR  zone  located  near  a 
loaded  crack  tip  of  ductile  materials.  Attempts  to  solve  for  the  stresses  and  strains 
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within  metals  exhibiting  power  law  hardening  proved  to  be  intractable  in  closed  form 
solution  and  finite  element  analysis  was  then  used  to  solve  these  types  of  problems. 

In  the  past,  the  path  independence  of  J  has  been  established,  excluding  contours 
that  pass  through  dissimilar  media,  holes  or  drastic  dimensional  changes  or  where  a 
large  amount  of  crack  extension  takes  place  [32].  The  stress  intensity  factors  K  and  J 
are  only  applicable  if  they  uniquely  characterize  o^/and  e,y  at  the  location  of  fracture 
initiation  [33],  It  has  been  shown  analytically  that  the  J-dominated  HRR  region  is  valid 
up  to  0.25  times  the  length  of  the  plastic  zone  radius  [34-37].  Thus  the  J-Integral  has 
been  used  to  study  fracture  problems  involving  plasticity  beyond  the  scope  of  LEFM 
based  on  the  K  concept  [38]. 

An  alternative  fracture  criterion  has  been  proposed  by  Sih  [39-40]  which  is 
capable  of  predicting  the  critical  crack  initiation  load  as  well  as  the  crack  propagation 
direction.  By  measuring  the  amount  of  strain  energy  absorbed  by  a  material  in  front  of  a 
crack  tip,  the  strain  energy  density  factors  Smin  and  Snua  may  be  calculated,  which  refer 
to  dilatation  (local  volume  changes)  and  distorticHi  respectively.  The  critical  value  of 
Smin  may  then  be  used  as  a  failure  criterion  to  prredict  the  failure  load  of  sharply  defined 
cracks  or  elliptical  notches. 

Toughness  versus  Work  of  Fracture  in 
Metal  Matrix  Composites 

Several  authors  [41-42]  have  pointed  out  that  characterizing  the  plane  strain 
fracture  toughness  (Kc  or  K/c)  in  composite  materials  is  not  as  straight  foru'ard  as  in 
isotropic  metals.  The  inherent  anisotropy  of  composite  materials  makes  the  fracture 
mode  mixed,  even  under  simple  tension.  ASTM  has  established  rigorous  testing 
procedures  (E399)  to  measure  the  inherent  fracture  toughness  of  quasi-brittle  materials; 
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however,  a  standardized  test  method  for  MMCs  has  not  yet  been  developed.  Friend 
[41]  points  out  that  Kc  is  not  a  measure  of  the  material's  homogeneous  fracture 
toughness  since  continuous  fiber  metal  matrix  composites  are  highly  anisotropic. 
Fracture  toughness  will  be  used  in  reference  to  the  parameter  rigOTOusly  established  by 
ASTM  E399  and  "toughness"  will  be  used  in  reference  to  a  measure  of  the  materials 
ability  to  absorb  energy. 

An  alternative  method  of  measuring  the  "toughness"  of  composite  materials  is  by 
measuring  the  fracture  energy  or  work  of  fracture.  The  following  definition  has  been 
given  by  Davidson  [42], ", . .  work  of  fracture--the  dissipation  of  the  elastic  energy 
imposed  in  the  specimen  by  the  externally  applied  load  per  unit  of  new  crack  surface 
formed."  This  measure  of  toughness  may  be  obtained  by  calculating  the  area  under  a 
load-crack  opening  displacement  (COD)  curve  during  3  point  bend,  or  the  area  under  the 
force-time  curve  in  an  impact  test.  Magata  and  Hall  [43]  have  measured  the  area  under 
the  load-CQD  curve,  and  found  radically  different  values  of  toughness  as  compared  to 
Kq  for  their  a-Al2C>3/Mg  system.  Specifically,  Kq  is  defined  as  the  experimentally 
determined  value  of  the  toughness  of  the  composite  that  does  not  satisfy  ASTM  standard 
E399.  For  composite  materials  it  may  be  appropriate  to  measure  the  toughness  using 
both  methods. 

Energy  Absorbing  Micromechanisms 

Pictured  in  Fig.  1. 1  is  the  crack  tip  process  zone  or  "local  heterogeneous  region" 
in  a  continuous  fiber  MMC.  Within  the  schematic  are  numerous  mechanisms  that 
contribute  to  the  total  fracture  in  a  fiber  composite  containing  a  distinct  interfacial  region. 
The  mechanisms  are  as  follows;  1)  fiber  pull-out,  2)  crack  bridging,  3)  crack  deflection, 
4)  fiber  breakage,  and  5)  matrix  yielding  and  plasticity.  Each  mechanism  will  be  treated 
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Breakage  Bridging 


Figure  1.1.  Schematic  of  the  crack  tip  process  zone  in  a  continuous  fiber  metal  matrix 
composite  showing  distinct  energy  absorbing  mechanisms. 


separately  in  detail,  with  the  exception  of  matrix  yielding  and  plasticity,  which  has  been 
previously  discussed. 
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Fiber  pull-out 

Additional  energy  is  expended  during  fiber  pull-out,  thus  by  increasing  the  pull¬ 
out  length,  additicnal  toughness  may  be  gained.  In  the  analysis  given  by  Taya  and 
Arsenault  [44],  the  average  pull-out  energy  per  unit  area,  is  given  as  follows: 
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The  maximum  pull-out  energy  is  given  when  /=/c*: 
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Equation  1.7  shows  that  the  maximum  energy  is  expended  in  pulling  out  fibers  given 
high  volume  fractions,  large  diameter  fibers,  high  fiber  strengths  and  low  interfacial 
strengths. 

Due  to  the  desire  to  retain  strength  and  stiffness,  the  bonding  between  fiber  and 
matrix  is  relatively  high.  Systems  such  as  a-Al203/Mg  [43],  Boron/Al  [45-46]  and 
SiC/Al  [47]  possess  good  chemical  bonding  and  experience  mechanical  clamping  due  to 
the  differences  in  thermal  expansion  coefficients  between  the  fiber  and  matrix.  These 
composites  are  not  likely  to  exhibit  fiber  pullout.  However,  composites  such  as  SCS- 
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6/Ti-15-3  exhibit  extremely  poor  interfacial  bonding  and  large  amounts  of  fiber  pull-out 
[48].  To  promote  fiber  pull-out,  Magata  and  Hall  [43]  thermally  cycled  their  a- 
Al20s/Mg  composite  to  debond  the  fibers  from  the  matrix  consequently  lowering  the 
degree  of  interfacial  bonding.  Although  the  measured  values  of  Kq  decreased  from 
22.7  to  9.1  MPaVm  when  tested  in  the  longitudinal  directicm,  the  work-of-fracture 
increased  from  0.202  to  0.354  J  from  the  uncycled  to  the  cycled  conditions.  The  75% 
increase  in  work-of-fracture  would  have  gone  unnoticed  if  both  experiments  had  not 
been  conducted. 

Crack  bridging 

In  the  case  when  the  matrix  separates  ahead  d*  the  crack  tip  and  the  fibers  are 
exposed  or  "bridge”  the  crack,  additicmal  energy  is  required  for  crack  propagation  [44]. 
The  additional  plastic  work  (Wp)  perftxmed  on  the  matrix  by  the  fibers  during  crack 
bridging  is  given  by: 

W^icd\af,-aj)b!n  <>-8) 

This  phenomenon  was  highly  noticeable  in  an  AI2Q3/AI  composite  formed  through  the 
Lanxide®  process  [49].  Although  the  reinforcement  was  discontinuous,  crack  bridging 
added  to  the  toughness  in  the  presence  of  lower  volume  fractions  of  alumina  In  order 
for  the  composite  to  exhilnt  crack  bridging,  a  favorable  stress  state  must  be  present 
Increasing  triaxiality  would,  in  effect  favor  the  rapid  nucleation  of  microvoids. 
Awerbuch  and  Hahn  [46]  noted  that  the  presence  of  bwon  fibers  in  an  aluminum  matrix 
acted  to  increase  the  triaxial  stress  state  ahead  of  the  crack  tip,  which  could  be  seen  frcHn 
the  large  regions  of  plastic  deformation. 
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Crack  deflection 

Irregular  crack  paths  expend  greater  energy  than  self-similar  cracks  due  to  the 
decreased  stress  concentration  at  the  crack  tip.  This  may  be  directly  caused  by  deflected 
crack  propagation  in  angle  ply  composites  or  interfacial  crack  propagation  in 
unidirectional  composites.  An  indirect  effect  erf  the  irregular  crack  path  is  crack  closure, 
which  shields  the  crack  tip  thus  reducing  the  effective  stress  concentration.  The 
irregular  crack  p>aths  significantly  increase  the  toughness  in  B/Al  [45]  and  a-Al2C^/Al 
[43].  Mahulikar  et  al.  [2]  have  shown  that  fiber  orientation  plays  an  important  role  on 
the  failure  mode  of  B4C  and  SiCyTi-6-4  metal  matrix  composites.  Self-similar  crack 
propagation  consisting  of  transfibral  and  interfacial  failures  occurs  when  the  crack  is 
oriented  at  an  angle  greater  than  45*  to  the  fiber  direction.  If  the  pre^rties  of  the 
interface  are  altered  by  thermal  cycling,  then  failure  occurs  primarily  along  the  interface. 
Thermal  cycling  forms  reaction  products  altwig  the  interface  creating  a  "band  <rf 
weakness"  in  which  failure  is  favored.  Clearly  the  degree  of  interfacial  bonding  is 
related  to  the  amount  of  crack  deflection. 

Recent  analytical  work  by  Gupta  and  Argon  [14]  cm  brittle  compxsite  interfaces 
suggests  that  fiber  breakage  or  delamination  can  be  controlled  by  varying  the  elastic 
constants  of  both  the  interphase  and  the  fiber.  This  work  has  been  made  possible 
through  new  techniques  of  measuring  the  material  propjerties  of  the  fiber  and  interphase 
using  a  unique  laser  spxdlation  technique  [50].  The  earliest  reference  to  this  typse  of 
work  can  be  traced  back  to  studies  of  cracks  adjacent  to  bimaterial  interfaces.  Gupta  and 
Argtm  use  the  following  pjarameters  A  and  p  to  define  the  degree  of  anisotropy  in  these 
composites: 


A  =  ^li 
S22 


(1.9) 
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These  parameters  (A  and  p)  form  the  basis  of  the  Dundur's  parameters  a  and  /3  [14]  as 
is  given  below: 
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The  Dundurs  parameters  p  and  p  rem^n  relatively  constant  when  values  of  the  elastic 
constants  for  typical  engineering  materials  are  substituted.  The  final  dependence  for  the 
generalized  delamination  chart  falls  onto  a  choosing  P7=l,  P2=l>  ^2=1  ^2=1.  By 

plotting  a  against  the  ratio  of  the  applied  transverse  stress  to  the  strength  of  the 
interface,  the  failure  mode  of  the  material  may  be  predicted  (see  Fig  1.2).  Above  the 
constructed  line,  cracks  will  deflect  along  the  interface,  while  the  lower  region  will  result 
in  fiber  failure.  Knowing  the  elastic  ccmstants  of  the  interface  and  the  fiber,  the  ratio  d* 
the  transverse  strength  to  the  longitudinal  strength  may  be  calculated.  In  this  manner, 
delamination  charts  were  analytically  formulated  to  engineer  composite  materials  that 
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provide  both  high  transverse  and  longitudinal  su-ength  as  well  as  high  toughness  (fiber 
pull-out).  Through  such  processes  as  CVD  and  PVD,  the  elastic  properties  of  the 
interface  may  be  tailored  to  suit  a  specific  application. 

Fiber  breakage 

The  ultimate  strength  of  a  CFM  is  highly  dependent  upon  the  strength  of  the 
individual  fibers.  The  strength  of  these  fibers,  unlike  that  of  the  metal  matrices,  is 
statistical  in  nature  and  may  be  described  using  Weibull  statistics.  The  less  scatter  in 
strength  the  fibers  possess,  the  easier  it  is  to  predict  the  failure  of  a  fibrous  composite. 
Such  a  distribution  can  be  given  by  the  following  relation  [44]; 

=  la^  exp{-laol]  (1.15) 


Micromechanical  Modeling 

Micromechanical  modeling  can  be  broken  into  two  areas  of  interest:  analytical 
modeling  and  finite  element  analysis.  Analytical  modeling  makes  simplifying 
assumptions  of  the  material's  behavior  yielding  constitutive  equations  in  closed  form. 
Finite  element  analysis  solves  for  the  stresses  and  strains  of  a  discretized  model 
numerically.  Both  methods  have  their  advantages  and  disadvantages.  Analytical  models 
provide  mathematically  simple  equations  that  can  be  used  for  design,  but  often  at  the 
expense  of  realism.  Under  careful  scrutiny,  finite  element  analysis  can  yield  reliable 
solutions,  but  may  require  large  amounts  of  computer  resources. 
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Analytical  Formulations 

To  date,  there  has  been  a  significant  amount  of  work  performed  on  determining 
the  mechanical  properties  of  composites  from  the  properties  of  the  constituents.  Hashin 
[51]  has  categorized  analytical  procedures  into  three  areas:  the  direct  approach,  the 
variational  approach  and  approximation  techniques.  These  procedures  may  be  applied  to 
statistically  isotropic  (particulate)  composites  as  well  as  to  fibrous  composites.  Of  these 
approaches,  only  the  models  using  the  direct  approach  will  be  considered.  The  direct 
approach  yields  an  exact  solution  of  the  effective  properties  of  a  composite  considering  a 
geometric  model.  Below  are  highlighted  some  popular  micromechanical  models  that  are 
relevant  to  metal  matrix  composites. 

Vanishing  fiber  diameter  model 

This  model  was  developed  by  Dvorak  and  Bahei-El-Din  [52,53]  specifically  for 
analyzing  micromechanical  problems  in  metal  matrix  composites.  The  model  was 
evaluated  using  the  2-dimensional  laminate  AGLPLY  program  by  Bigelow  and  Naik 
[15].  The  values  of  the  longitudinal  modulus  and  Poisson's  ratio  were  found  to  be 
within  6%  of  the  experimentally  determined  values  in  the  B/Al  and  SiC/Ti  [0]8 
laminates.  Although  AGLPLY  was  designed  to  predict  the  laminate  response, 
individual  constituent  stresses  can  be  calculated.  It  was  found  that  the  constituent 
stresses  developed  by  mechanical  loads  were  more  accurate  than  stresses  developed  by 
thermal  loads  when  compared  to  finite  element  calculations  [15]. 

Aboudi  model 

The  Aboudi  model  [54,55]  is  capable  of  describing  the  effective  behavior  of 
elastoplastic  composites.  The  periodic  and  geometrical  arrangements  of  rectangular 
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parallelepipeds  in  the  Aboudi  model  allow  both  continuous  fiber  and  particulate 
composites  to  be  studied.  The  nonlinearity  of  the  metal  matrix  has  been  accounted  for 
using  elastoplastic  temperature  dependent  constitutive  relations.  The  complex  derivation 
and  results  are  too  lengthy  to  discuss  in  detail  here.  For  the  complete  deri\'ation  and 
constitutive  relations,  see  reference  54. 


Cox  model 

The  Cox  model,  as  summarized  by  Johnson  and  Birt  '56],  was  designed  to 
analyze  short  fiber  systems.  This  model  ac  xjunts  for  the  shear  and  axial  stress 
distribution  along  an  embedded  *!t  jr.  The  Cox  model  also  derives  the  longitudinal 
stiffness  as  given  in  equations  (1.16-1.20)  according  to  the  coordinate  system  shown  in 
Fig.  1.3. 
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In  their  evaluation,  Johnson  and  Birt  found  the  Cox  model  accurately  predicted 
longitudinal  stiffness.  However,  the  Cox  model  should  not  be  used  for  very  short  fiber 
composites  (1/d^)  or  for  particulate  composites  due  to  inadequate  load  transfer  through 
shear. 


Figure  1.3.  Coordinate  system  used  in  deriving  longitudinal  stiffness  of  the  Cox 
model. 


Paul  model 

The  Paul  model  was  derived  to  predict  the  elastic  modulus  of  particulate 
composite  materials.  This  model  assumes  that  the  Poisson's  ratio  of  the  matrix  equals 
that  of  the  reinforcement  and  that  the  strain  in  the  matrix  equals  the  strain  seen  by  the 
particles.  The  transverse  and  longitudinal  elastic  moduli  are  given  in  equation  (1.21). 
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Johnson  and  Bin  [56]  found  that  the  Paul  model  predicted  an  upper  bound  of  elastic 
moduli  in  particulate  metal  matrix  composites,  but  should  not  be  used  for  whisker 
reinforced  metals. 


Computational  Finite  Element  Analysis 

The  finite  element  method  (FEM)  has  been  the  solution  procedure  of  choice  for 
complex  problems  involving  mechanical  and  hygrothermal  loading  of  composite 
materials  [48,57-64].  This  method  has  been  applied  to  brittle  [57,58]  and  ductile  matrix 
[59-63]  composites.  Current  and  past  computational  studies  of  composite  materials  can 
be  separated  into  two  levels:  a  continuum  or  global  level  and  a  micromechanical  or  local 
level.  Some  investigators  have  chosen  to  link  these  regimes  through  a  local-global 
analysis  where  the  global  stress  state  and  loading  history  dictates  the  boundary 
conditions  for  the  local  model.  In  the  following  sections,  both  levels  of  analysis  will  be 
described  along  with  example  case  studies. 

Continuum  approach 

The  object  of  a  global  analysis  is  to  describe  the  overall  mechanical  response  of  a 
composite  structure.  The  material  is  considered  statistically  homogeneous  down  to  the 
ply  level.  Several  finite  element  codes  have  been  written  to  account  for  fibrous 
reinforcement  and  laminate  design. 
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Aboudi's  vanishing  fiber  diameter  (VFD)  model  and  the  plasticity  theories 
dev  eloped  by  Dvorak  and  Bahei-El-Din  [52,53]  inspired  codes  such  as  PAC78  and 
PAFAC.  NASA  researchers  [6,15,62]  have  used  these  codes  to  observe  strain 
distributions  ahead  of  geometrical  discontinuities  (such  as  notches  and  holes). 
Relatively  good  comparisons  with  experimental  results  were  noted. 

The  three-dimensional  package,  PAFAC,  has  been  extensively  used  in  studying 
boron/aluminum  and  silicon  carbide/aluminum  systems  [15,62].  The  eight  noded 
hexahedral  elements  represent  a  composite  continuum  incorporating  matrix  and  fiber 
failure  theories.  The  following  equation  predicts  fiber  failure  in  PAFAC; 
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Where  ai  i  is  the  tensile  stress  in  the  fiber  direction  and  012  is  the  shear  stress. 

An  ultimate  strain  criterion  is  used  to  calculate  matrix  failure: 

(^f + (^nf  -  + 3(.r2N  (1.23) 

In  a  study  of  SCS-2/6061  Al,  Johnson  and  Bigelow  [62]  noted  good  prediction  of  the 
normalized  strains  ahead  of  a  slit  in  [0/90]2s  laminates.  The  numerical  predictions  from 
PAFAC  were  compared  to  results  obtained  with  an  array  of  five  strain  gages  positioned 
ahead  of  the  notch.  In  [±4^2s  laminate  poor  correlation  was  found,  possibly  due  to 
fiber  rotation  and  interlaminar  yielding. 

To  overcome  possible  interlaminar  nonlinearity,  Hinrichsen  and  Palozotto  [65] 


have  devised  a  cubic  spline  function.  In  effect,  the  assumption  of  plane  sections 
remaining  plane  is  eliminated.  This  has  been  a  problem  in  thick  plate  applications.  In 
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the  past,  thin  shell  elements  ha\'e  been  "artificially  softened"  to  simulate  real  life 
applications,  whereas  Hinrichsen  and  Palozotto  have  generated  accurate  results  without 
empirically  reducing  G13  and  G23. 

A  final  application  of  continuum  FEA  concerns  predicting  the  mechanical 
properties  of  laminates.  Johnson  et  al.  [48]  have  used  AGLPLY,  a  code  employing 
Aboudi's  plasticity  theories,  for  predicting  the  elastic/plastic  behavior  of  symmetric 
c(»nposite  plates  supporting  in-plane  loads.  AGLPLY  has  successfully  predicted  overall 
laminate  moduli,  ply  stresses  and  strains  throughout  the  loading  history. 

Micromechanics  approach 

A  mcM-e  desirable  approach  would  be  to  predict  the  composites  mechanical 
properties  based  only  on  the  mechanical  properties  of  the  individual  constituents.  In  the 
next  segment,  three  representative  micromechanics  approaches  will  be  discussed.  Each 
uses  the  material  properties  of  distinct  fiber  and  matrix  regions.  Dvorak  and  Teply  [61] 
developed  an  ingenious  energy  balance  method  to  predict  instantaneous  composite 
stiffness,  while  Lerch  and  Melis  [66]  used  a  straight  ahead  approach  to  describe  the 
stress  state  and  slip  band  regions  in  cross  ply  laminates.  Pagano  and  Tandon  [64,67] 
utilized  a  3-cylinder  assemblage  model  to  study  the  effects  erf  various  coatings  on  the 
stress  distribution  around  a  thermally  loaded  fiber. 

On  a  micromechanical  scale,  it  is  important  to  establish  a  working  representative 
volume  element  (RVE)  which  accurately  represents  the  composite.  The  only  models  for 
which  mathematically  precise  solutions  may  be  obtained  are  the  composite  cylinder 
assemblage  and  periodic  anays  of  identical  fibers.  The  former  may  be  solved 
analytically,  while  the  latter  must  be  solved  numerically.  To  date,  a  variety  of  periodic 
fiber  arrays  have  been  used  to  determine  mechanical  property  data  for  composite 
materials. 
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Dvorak  and  Teply  [61]  have  presented  a  guideline  that  should  be  used  to 
construct  accurate  representative  volume  elements:  a)  when  repeated  it  continuously 
covers  the  entire  macroscopic  volume  of  the  composite  and  b)  when  loaded  by  uniform 
stresses  or  strains,  the  local  stresses  and  strains  are  identical  within  each  RVE  in  the 
macroscopic  volume.  Likewise,  they  represent  the  macroscopic  composite  by  an 
equivalent  homogeneous  volume  (EHV)  element.  By  equating  the  energy  changes  in  the 
RVE  and  the  EHV  when  subjected  to  the  identical  boundary  conditions,  they  determine 
the  instantaneous  stiffness  of  the  EHV.  The  energy  change  (AIT)  from  a  known 
solution  state  to  a  slightly  perturbed  state  has  been  derived  for  the  RVE  and  the  EHV: 

AHemv  =  An  Aa^B^LBAa  -  Ao^BAa  (1.24) 

2 

AHrve  =  An  =  I  Aa^KAa  -  A^^Aa  (1.25) 

2 

Through  boundary  cemditions  applied  to  the  RVE  and  EHV,  Dvorak  and  Teply 
accurately  predicted  upper  and  lower  bounds  of  elastic  moduli  using  ABAQUS. 
Displacement  or  strain  was  applied  to  the  EHV  to  give  an  upper  bound  estimate,  while 
stress  or  traction  was  tqjplied  to  give  a  lower  bound  estimate.  The  instantaneous 
stiffness  L  of  the  macroscopic  composite  is  at  this  point  unknow  n.  Next,  the  identical 
boundary  conditions  were  applied  to  the  RVE  having  know  n  instantaneous  stiffness  L, 

(i  represents  the  fiber  and  matrix).  The  energy  terms  in  both  the  local  and  global  models 
were  then  equated  to  determine  Lgiobal-  The  following  formula  was  used  by  Dvorak 
and  Teply  to  determine  the  energ)'  change,  ATI,  for  an  upper  bound  estimate  using  a 
displacement  approximation: 
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Ati  = 


tdudS 


(1.26) 


The  equilibrium  model  was  then  used  for  determining  a  lower  bound  estimate: 


(1.27) 


Energy  principles  were  used  to  show  the  u^^r  and  lower  bound  estimates.  Dvorak  and 
Te|4y  utilized  the  maximum  principle  for  strain  rates  on  the  upper  bound  while  the 
elastic  minimum  ccmi|rfemcntary  energy  theorem  was  used  for  the  lower  bound. 

Lerch  and  Melis  [66]  predicted  the  stress  strain  response  of  cross-ply  laminates 
from  a  micromechanical  model.  Residual  stresses  and  fiber  debonding  were  accounted 
fcM*  through  gap  elements  that  opened  at  a  critical,  predetermined  stress  state.  The 
titanium  matrix  was  modeled  as  elastic/plastic,  while  the  fibers  were  assumed  to  remain 
elastic  to  failure.  Using  the  MARC  finite  element  package,  good  correlation  of  the  axial 
and  transverse  stress  strain  behavior  was  found  between  the  model  and  experimental 
results. 

A  three  phase  cylinder  model,  distinguishing  fiber,  matrix  and  interphase  regitxis 
was  utilized  by  Pagano  and  Tandon  [67]  to  study  the  elastic  response  of  multi- 
directitmal  coated  liber  composites.  By  parameterizing  the  mechanical  (xt^rties  and 
geonetry  of  the  interphase,  the  stress  intensity  about  a  thermally  excited  fiber  was 
documented  experimentally  and  numerically.  Their  wcM’k  centered  on  Nicalon  reinforced 


BMAS,  with  fiber  coatings  composed  of  nickel  and  carbon.  The  derived  constitutive 
relations  were  solved  using  the  numerical  package,  NDSANDS. 
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Objective  and  Scope  of  the  Investigation 

The  primary  objective  of  this  investigation  was  to  develop  a  local-global 
micromechanical  model  for  continuous  fiber  metal  matrix  composites.  This  finite 
element  model  is  to  account  for  energy  absorbing  mechanisms  such  as  fiber  debonding 
and  matrix  plasticity  on  a  local  level.  Parametric  studies  of  the  interfacial  bond  strength 
and  fiber  vdume  fraction  are  to  be  presented  and  analyzed  to  provide  insight  on  the 
mechanics  of  composite  fracture.  It  was  also  intended  that  fracture  in  monolithic  metals 
should  be  well  understood  and  could  be  related  to  matrix  failure  in  the  composite  model. 

Chapter  2  will  describe  the  finite  eluent  analysis  procedure  and  will  outline  the 
various  case  studies.  Modifications  that  were  made  to  the  PAFST  code  will  be 
highlighted  along  with  various  post  processors  used  to  analyze  the  data.  Chapter  3  will 
discuss  the  oomputational  results  of  both  the  two  dimensional  elastic-plastic  fracture 
proUems  as  well  as  the  local-global  composite  fracture  problem.  Chapters  4  and  5  will 
present  impcmant  conclusions  and  give  recommendations  for  future  work,  respectively. 


CHAPTER  2 


NUMERICAL  PROCEDURE 


To  study  the  fracture  characteristics  of  metal  matrix  composites,  a  two  part 
investigation  was  necessary.  First  a  detailed,  nonlinear  finite  element  analysis  has  been 
carried  out  on  a  variety  of  structural  alloys.  The  numerical  code,  PAPST  [68],  was 
modified  to  peiform  the  present  analysis.  Furthermore,  a  local-global  investigation  has 
been  conducted  using  the  commercially  available  package,  ANSYS. 


Fracture  of  Monolithic  Metals 

Listed  below  are  the  aq»tHlities  and  modifications  that  have  been  made  to  the 
PAPST  code.  The  execution  procedure  and  additional  postprocessors  have  been 
documented  to  account  for  all  the  changes. 

Capabilities  of  PAPST 

PAPST  (Plastic  Axisymmetric/Planar  Structures)  can  be  traced  back  to  another 
finite  element  code,  APES,  which  was  limited  to  linear  clastic  fracture  problems 
[69,70].  Gifford  and  Nash  subsequently  modified  this  code  to  accept  nonlinear  material 
characteristics  and  have  published  their  findings  on  crack  problems  in  elastic-plastic 
materials  [71].  To  handle  the  mathematical  instability  at  the  crack  tip,  a  special  enriched 
mode  I  singularity  element  was  incorporated  with  8x8  Gaussian  quadrature. 
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Away  from  the  highly  stressed  regions,  PAPST  uses  12-noded  isoparametric 
quadrilateral  2-D  elements  with  4x4  Gaussian  quadrature.  All  of  the  nodc^  are  located 
on  the  perimeter  of  the  element.  These  elements  allow  stresses  to  vary  over  the  face  of 
the  element  thus  requiring  fewer  elements  to  model  a  given  structure.  Both  element 
types  are  capable  of  handling  linear  hardening  materials  such  as  those  studied  by 
Hutchinson  [30],  Rice  and  Rosengren  [31].  A  comparison  to  constant  stress  (strain) 
triangular  elements  is  given  in  the  APES  manual  [69,70]  showing  the  accuracy  of  using 
these  enhanced  elements. 

Solutions  are  obtained  nonlinearly  using  the  Newton-Raphson  iteration 
technique  [72].  The  nonlinear  capability  of  PAPST  incorporates  the  J2-incremental  flow 
thetxy  of  plasticity  and  the  Von  Mises  yield  criteria.  The  incremental  plasticity  theory  is 
given  as  [68]: 


Cij  = 


1  + V 


fo’’ 


(2.1) 


i  -  1  +  V  ; 

^•J - £ — 

(otherwise)  (2.2) 

where  is  the  effective  stress  and  Oe  is  the  Von  Mises  effective  stress  with  respect  to 
the  coordiirates  of  the  current  center  of  the  yield  surface.  PAPST  approximates  the  true 
stress-strain  curve  of  a  material  using  either  multilinear  or  Ramberg-Osgood  power 
hardening  approximations.  PAPST  is  capable  of  modeling  either  purely  kinematic  or 
isotropic  hardening  materials  or  may  model  a  mixture  of  the  two.  Kinematic  hardening 
is  when  the  center  of  the  yield  surface  is  found  to  move,  while  isotropic  hardening  is 
when  the  yield  surface  expands  uniformly.  PAPST  handles  both  plane  stress  and  plane 


strain  conditions.  TheJ-Integral,  an  important  parameter  in  nonlinear  fracture 
mechanics,  may  be  calculated  over  10  paths  surrounding  the  crack  tip. 


Modifications  Made  to  PAPST 

The  PAPST  code  has  been  in  existence  in  some  form  or  another  since  1972.  In 
1984  this  code  was  modified  to  run  on  the  VAX  mainframe  environment  by  Bhagat  et 
al.  [73].  The  code  was  subsequently  termed  ARLPAPST  which  was  substantially 
expanded  in  its  solution  capability  and  integrated  pre-  and  postprocessing  subroutines. 
ARLPAPST  handled,  at  most,  50  elements  or  500  nodes  as  was  the  case  with  PAPST. 
The  original  inLitt  of  this  thesis  was  to  model  composite  materials  using  ARLPAPST, 
therefore  the  code  was  altered  to  handle  an  order  of  magnitude  greater  number  of 
elements  and  nodes. 

Several  arrays  including  those  controlling  boundary  conditicms,  thermal  effects, 
spring  elements  and  global  solution  matrices  were  altered  in  size.  From  personal 
communications  with  Gifford  [68],  the  following  formula  w'as  used  to  determine  the 
array  size  of  G,  the  global  solution  matrix,  where  N  is  the  number  of  nodes: 

G  =  2*(N+1,^  (2.3) 

Additional  Postprocessors 

Modifications  were  made  to  output  the  strain  energy  density,  selected  7-IntegraI 
values  and  plots  of  the  deformed  mesh.  To  analyze  crack  instability,  values  of  the  strain 
energy  density  (dW/dV)  were  required.  Since  ARLPAPST  generates  dWIdV  w'hen 
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calculating  J,  a  write  statement  placed  in  the  ccxle  outputs  dW/dV  w  ithout  additional 
effort.  This  write  statement  w'as  designed  to  only  print  out  dW/dV  along  the  first  J- 
Integral  path  to  save  disk  space.  The  nodes  along  the  path  do  not  need  to  form  a  closed 
contour.  The  post  processors  that  were  written  to  cut  J  and  generate  postscript  output  of 
the  deformed  structure  are  listed  in  appendices  A  and  B  respectively.  Finally,  contour 
plots  were  generated  using  SAS  graph  on  IBM  3179  terminals  and  pen  plotters. 

Solution  Procedure 


The  solution  procedure  for  analyzing  nonlinear  elastoplastic  fracture  problems 
with  ARLPAPST  is  outlined  in  Fig.  2.1. 

Materials  Under  Investigation 

Two  types  of  steel  (HY-130  and  HY-140)  as  well  as  two  types  of  aluminum 
alloys  (7075  and  5083-0)  have  been  investigated.  The  compact  tension  specimens  that 
have  been  analyzed  conform  to  current  ASTM  standards,  see  Table  2.1.  Multilinear 
material  modeling  of  the  true  stress  strain  curves  are  used  to  represent  all  the  materials  as 
is  shown  in  Fig.  2.2.  Typical  mesh  designs  are  presented  in  Chapter  3  (refer  to  Figs. 
3.1  and  3.2). 


Local-Global  Micromechanical  Modeling 


The  second  stage  of  this  investigation  centered  on  the  fracture  behavior  of  SCS- 
6yTi-15-3  composite.  Johnson  et  al.  [48]  performed  exhaustive  experimental  work  on 
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Figure  2. 1.  Flowchart  for  analyzing  elastic-plastic  fracture  problems  using  the 
modified  ARLPAPST  finite  element  code 
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Table  2.1.  Compact  tension  specimen  geometry  and  material  data 


Material 

IVCmm) 

fl(mm) 

a(mm) 

aJW 

(MPa) 

Condition 

Reference 

HY-130 

50.8 

25.4 

39.04 

0.7685 

896.3 

plane  e 

[71] 

50.8 

25.4 

41.22 

0.8115 

50.8 

25.4 

42.92 

0.8843 

50.8 

25.4 

25.40 

0.5000 

HY-140 

50.8 

25.4 

39.04 

0.7685 

965.3 

plane  e 

[711 

7075  A1 

50.8 

25.4 

39.04 

0.7685 

413.7 

plane  e 

[74] 

50.8 

25,4 

41.22 

0.8115 

5083-0  A1 

152.4 

45.7 

81.28 

0.5333 

121.8 

plane  a 

[75] 

152.4 

45.7 

81.28 

0.5333 

121.9 

plane  e 

True  stress-strain  cun’es  of  HY-130  steel,  HY-140  steel  and  5083-0 
aluminum  and  7075  aluminum. 
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this  material,  characterizing  the  stress  strain  behavior  of  multi-ply  systems.  Their  work 
indicates  that  SCS-6/Ti-15-3  possesses  a  weak  interface  (138  MPa)  which  they  have 
documented.  Lerch  and  Melis  [66]  have  used  their  work  to  mode!  cross-ply  composites 
and  consider  the  work  of  Johnson  et  al.  [48]  to  be  comprehensive. 

Material  Properties 

The  compiled  properties  of  silicon  carbide  SCS-6  fibers,  as  fabricated  Ti-15V- 
3Cr-3Al-3Sn  and  unidirectional  as  fabricated  SCS-6/Ti-15-3  composite  are  listed  in 
Table  2.2.  The  nonlinear  stress-strain  data  used  for  the  local  model  is  presented  in  Fig. 
2.3.  Only  the  elastic  properties  of  the  composite  are  used  on  the  global  level. 

Numerical  Approach 

First  a  global  elastic  analysis  was  performed  on  composite  using  specimen 
dimensions  extracted  ftx)m  the  work  of  Johnson  et  al.  [48].  The  center  cracked  plate 
under  study  was  modeled  using  the  boundary  conditions  shown  in  Fig.  2.4.  A  fully 
three  dimensional  global  analysis  was  not  performed  due  to  a  lack  of  experimental  data 
of  the  necessary  five  independent  elastic  constants.  Instead  a  two-dimensional,  plane 
strain  model,  which  incorporates  three  independent  elastic  constants,  was  studied  with 
orthotropic  elements  (STIF  82)  using  the  ANSYS  finite  element  code.  A  uniform 
displacement  of  2.5  pm  was  applied  to  the  top  edge  of  the  model,  corresponding  to  the 
minimum  displacement  necessary  to  induce  a  transverse  stress  greater  than  138  MPa,  the 
stress  at  which  fiber  debonding  occurs.  Four  finite  element  meshes  with  varying 
number  of  elements  were  then  anah’zed  to  optimize  the  numerical  solution.  A  sample 
input  file  is  presented  in  Appendix  C.  The  results  of  this  optimization  will  be  presented 
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Table  2.2.  Compiled  properties  of  fiber,  matrix  and  composite  for  SCS-6/Ti-15-3 
system  (properties  extracted  from  reference  48). 


Properties  of  the  fiber 

Coefficient  of  Thermal 
Expansion  (a) 

2.70  X  10-6  T-i 

4.09  X  10-6  ‘c-l 

Poisson's  Ratio  (v) 

0.25 

0.25 

Fiber  Diameter  (d) 

0.0056  in 

0.014224  cm 

Longitudinal  Elastic 

Modulus 

58.0  X  106psi 

399.9  GPa 

Transverse  Elastic  Modulus 

58.0  X  106  psi 

399.9  GPa 

Ultimate  Strength 

465  X  103  psi 

3208  MPa 

Properties  of  the  matrix  (as  fabricated) 

Coefficient  of  Thermal 
Expansion  (a) 

5.4  X  10^  *F'l 

9.72  X  10-6  ‘c-l 

Poisson's  Ratio  (v) 

0,351 

0.351 

Longitudinal  Elastic 

Modulus 

13.4  X  106  psi 

92.39  GPa 

Yield  Stress 

100.0  X  ICP  psi 

689.5  MPa 

Ultimate  Stress 

135.4  X  10^  psi 

933.55  MPa 

Properties  of  the  composite  (as  fabricated  [Olg  laminate) 

Fiber  Volume  Fraction 

0.325 

0.325 

Poisson's  Ratio  (v) 

0.299 

0.299 

Debond  Stress 

20  X  10^  psi 

138  MPa 

220  X  10^  psi 


Ultimate  Strength 


1518  MPa 


3.175  mm 


3.175  mm 


Figure  2.4.  Schematic  of  the  macroscopic  model  showing  boundary  conditions,  fiber 
orientation  and  specimen  dimensions. 
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in  Chapter  3.  The  displacements  in  the  most  highly  stressed  region  of  the  optimized 
model  were  then  used  as  the  boundary  conditions  for  the  local  model. 

The  three-dimensional  local  model  discretizes  the  fiber,  matrix  and  intetphase 
into  a  unit  cell  after  the  work  of  Dvorak  and  Teply  [61].  A  representative  unit  cell  with  a 
volume  fraction  of  0.325  is  shown  in  Rg.  2.5.  Equation  2.4  enaUes  the  calculation  (rf 
the  interfiber  spacing  (X)  of  the  unit  cell  model. 


A  =  2r 


(2.4) 


The  fiber  and  matrix  regions  were  constructed  with  20-noded  isoparametric  elements 
(STIF  95)  using  the  ANSYS  finite  element  code,  while  gap  elements  were  placed  alcmg 
the  fiber-matrix  interface. 

A  parametric  study  of  the  interface  strength  was  performed  by  varying  the  fiber- 
matrix  debond  stress.  The  tensile  strength  of  the  interface,  c^,  was  set  equal  to  138 
MPa,  corresponding  to  the  experimentally  determined  debond  stress  [48].  As  an 
extreme  case,  oj  was  equated  to  <y,  the  fiber  strength.  From  the  initial  analysis  this  gave 
unreasonable  results  and  an  intermediate  value  was  chosen,  corresponding  to  the 
strength  of  the  matrix.  The  stiffness  of  the  gap  element  is  assumed  to  be  equal  to  that  of 
the  matrix  (92.4  MPa).  Various  interface  strengths  are  inccxporated  through  changes  in 
radial  interference  (with  units  of  displacement).  To  calculate  the  radial  strain  seen  by  the 
gap  elements  and  the  surrounding  material,  the  effective  length  acting  (xi  an  individual 
element  is  assumed  to  be  of  length  0.5A  (see  equation  2.4).  Thus  the  interface  strength 
is  modeled  as  an  effective  strain  acting  upon  the  fiber-matrix  interface.  Another  effective 
length  may  or  may  not  be  appropriate,  but  this  length  seems  to  be  reasonable.  Using 
Hooke's  law,  the  desired  interface  strength  is  divided  by  the  gap  stiffness  to  obtain  the 


237.6  ]*m 


71.12  /<m 


71.12  }tm 


0.01  /<m  (typ.) 


Figure  2.5.  Representative  unit  cell  distinguishing  fiber,  matrix  and  interphase  of  a 
0325  volume  fraction  SiC/Ti  composite.  The  fiber  is  modeled  as  elastic, 
while  the  matrix  is  elastoplastic.  Gap  elements  are  used  to  simulate  the 
interfacial  region. 
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gap  strain.  These  strain  values  are  multiplied  by  the  factor  0.5A  to  obtain  the 
interference  length.  For  model  continuity,  a  small  mathematical  separation  (0.01  pm)  is 
placed  between  the  tu'o  nodes  of  the  gap  element. 

The  pictorial  shown  in  Fig  2.6  emi^iasizes  how  the  boundary  conditions  of  the 
global  model  are  transferred  to  the  local  model.  Part  A  of  Fig.  2.6  represents  the  entire 
specimen  modeled  by  Jdmscm  et  al.  [48].  Through  the  use  of  symmetry'  conditions, 
part  B  of  Fig.  2.6  was  used  as  the  global  or  "macro*  model.  The  horizontal  and  vertical 
displacements  were  extracted  from  the  nodes  of  the  global  model  correspcmding  to  the 
tc^  and  bottom  layers  the  local  model  as  shown  in  part  C  of  Fig.  2.6.  The  height  of 
the  microscopic  model  (380.0  pm)  corresponds  to  approximately  four  layers  d[  "macro¬ 
model*  nodes  as  shown  in  Fig.  2.7.  These  displacements  were  plotted  and  fit  with  a 
pdyntxnial  curve.  The  fitted  curves  and  corresponding  dis(dacements  will  be  jxesented 
in  the  following  chapter.  When  the  correlatioi  of  the  data  points  was  less  than  0.92, 
then  visual  interpolation  from  the  graj^  was  used.  The  fitted  displacement  boundary 
conditions  were  then  transferred  to  the  3-dimensiona]  local  model  via  coOTdinate 
transformation  as  is  shown  in  part  D  of  Fig.  2.6.  The  (x',  y',  z')  axes  denote  the 
coordinate  system  used  to  construct  the  local  model.  The  (X,  Y,  Z)  axes  denote  the 
cocMdinate  system  of  the  global  model.  The  displacements  along  the  global  Z  directicm 
are  assumed  to  remain  unchanged. 
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Global  Level 


Bgure  2.6.  PictcMial  flowchart  showing  the  global  model  in  parts  A  and  B  and  the 
local  model  in  parts  C  and  D.  Note  that  the  global  model  is  2- 
dimensional  and  that  the  local  model  is  3-dimensional. 
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Figure  2.7,  Nodal  plot  of  the  global  model  near  the  crack  tip  region.  The  outlined 
section  represents  the  dimensions  of  the  local  model  superimposed  upon 
the  Tiguie. 
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CHAPTER  3 


RESULTS  AND  DISCUSSION 

The  objective  of  our  research  is  divided  into  tw'o  separate  stages.  First,  a 
detailed  analysis  of  elastic-plastic  fracture  in  unreinforced  metals  is  conducted,  followed 
by  a  study  of  the  micromechanics  of  fracture  in  silicon  carbide  reinforced  titanium.  The 
study  of  moncdithic  steel  and  aluminum  alloys  gave  the  opportunity  to  verify  work  that 
has  been  performed  in  the  past  as  well  as  serve  as  a  benchmark  for  problems  solved 
using  the  updated  version  of  ARLPAPST. 

Monolithic  Materiais 

Due  to  the  ductile  nature  of  steel  and  aluminum  we  were  interested  in  the  crack 
tip  plastic  zone  (CTPZ)  and  its  dependence  on  the  crack  length.  7-Integral  values  were 
established  on  several  different  paths  in  both  steel  and  aluminum  specimens  and  the 
relative  path  independence  of  J  will  be  defended.  The  strain  energy  parameters 
(dWldV)c  and  {dWldV),„i„  presented  by  Bhagat  et  al.  [72]  w'ere  verified  and  w  ill  be 
discussed  at  length.  Finally,  stress  and  strain  contour  plots  are  presented  to  study  how 
the  stresses  are  distributed  in  front  of  the  crack  tip. 
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Mesh  Deformation 


Plots  of  the  deformed  meshes  of  5083-0  aluminum  and  7075  aluminum  are 
given  in  Figs.  3.1  and  3.2,  respectively.  The  deformed  mesh  shown  in  Fig.  3.1 
corresponds  to  an  applied  load  of  6122.5  pounds,  while  that  in  Fig.  3.2  corresponds  to 
an  applied  load  of  2844  pounds.  The  J-lntegral  paths  are  located  as  shown.  These 
loads  represent  the  crack  initiation  loads.  In  the  next  section,  the  crack  initiation  load 
will  be  described  in  terms  of  the  strain  energy  density. 

Strain  Energy  Density 

The  strain  energy  density  (dW/dV)  in  front  of  the  crack  tip  as  described  by  Sih 
[39,40]  was  calculated  using  ARLPAPST.  At  a  critical  radius  r^,  the  strain  energy 
density  dropped  to  a  minimum  value  designated  (dWldV)„i„  as  is  shown  in  Figs.  3.3 
and  3.4.  Earlier  studies  have  shown  that  this  "U-shaped"  signature  is  found  only  in  the 
^=0*  direction  ahead  of  the  crack  tip  [72).  Note  that  in  Fig.  3.4  the  strain  energy 
density  exhibits  a  jagged  signature.  This  corresponds  to  unloading  at  the  crack  tip. 
These  values  of  (dWldV)mm  were  plotted  against  the  applied  load  and  at  a  critical  load, 
/**,  the  values  of  dWIdV  began  to  decrease  as  is  shown  in  Fig.  3.5.  This  peak  value  of 
the  strain  energy  density  signifies  the  point  at  which  the  structure  begins  to  expend 
energy  in  the  form  of  crack  growth.  The  maximum  value  of  (dW/dV)min  is  denoted  as 
(dW/dV)*,nin- 

Three  different  crack  lengths  in  HY-130  steel  and  two  crack  lengths  in  7075 
aluminum  were  analyzed.  For  the  specimens  in  which  the  numerical  solutions  remained 
stable,  the  failure  load  can  be  approximated  by  extrapolating  the  curve  down  the  load 
axis.  For  example,  in  HY-130  steel  with  o=1.623,  the  extrapolated  curve 
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Crack  tip 


Figure  3. 1 .  Picture  showing  the  deformed  upper  half  of  a  compact  tension  specimen 
model  constructed  of  5083-0  aluminum.  The  crack  tip  (a=3.199  in)  is 
located  as  shown. 


Figure  3.2.  Picture  showing  the  deformed  upper  half  of  a  compact  tension  specimen 
model  constructed  of  7075  aluminum.  The  crack  tip  (a= 1.537  in)  is 
located  as  shown. 
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Figure  3.3.  Strain  energy  density  plotted  against  the  radial  distance  from  the  crack  tip 
in  HY-130  steel  (0=1.537  in). 
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Figure  3.4.  Strain  energy  density  plotted  against  the  radial  distance  from  the  crack  tip 
in  5083  aluminum  (<7=3.199  in). 
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reaches  the  load  axis  at  4200  pounds,  and  when  a=  1.7623,  the  predicted  failure  load  is 
approximately  2100  pounds.  The  fracture  loads  of  the  aluminum  specimens  are  less  and 
the  trend  is  similar  to  that  observed  in  the  steel.  The  fracture  loads  of  the  smaller  crack 
lengths  are  consistent  with  fracture  mechanics  concepts. 

Crack  Tip  Plasticity 

The  principal  Von  Mises  stress  distribution  in  loaded  compact  tension  specimens 
of  5083  Al,  7075  A1  and  HY-130  steel  are  shown  in  Figs.  3.6, 3.7a,  3.7b,  3.8a,  3.8b, 
and  3.8c.  The  increasing  size  of  the  crack  tip  plastic  zone  with  decreasing 
crack  length  agrees  with  accepted  plasticity  theories.  For  example,  in  Figs  3.7a  and 
3.7b,  the  ligament  ahead  of  the  crack  tip  for  the  larger  crack  (a= 1.623)  exhibits  full 
plasticity  for  all  the  nodes  along  the  ligament,  while  for  the  smaller  crack  (a=1.537) 
there  is  a  small  region  along  the  ligament  that  remains  in  the  elastic  region. 

Even  though  all  the  J  contour  paths  in  7075  aluminum  (as  shown  in  Fig.  3.7a) 
pass  through  the  CTPZ  at  F*,  the  variation  from  the  averaged  value  of  J  is  only  3%. 
However,  when  the  load  is  increased  to  1. 13P*  the  variation  from  the  average  J  value 
increases  to  7%.  Figure  3.9  illustrates  how  J  varies  along  the  crack  ligament  in  7075 
aluminum  (<1=1.537)  for  load  increments  from  0.21P*  to  1.13P*.  Similar  trends  can  be 
observed  in  HY-130  and  7075  Al  as  shown  in  Figs.  3.10, 3.1 1  and  3.12. 

The  crack  initiation  load  may  be  predicted  by  studying  the  variation  of  J  as  a 
function  of  applied  load.  Typically  these  plots  exhibit  asymptotic  behavior,  approaching 
a  critical  value  of  P.  This  information  can  then  be  used  to  predict  the  critical  loading 
state  of  a  cracked  structure.  The  work  of  Hilton  and  Gifford  [71]  has  been  confirmed  in 
this  investigation  through  the  plot  shown  in  Fig  3. 13. 


Figure  3.6.  Contour  plot  of  the  Von  Mises  effective  stress  (psi)  in  5083-0  aluminum 
(a=3. 199  in).  The  model  is  shown  to  scale  in  inches. 
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Figure  3.8a. 


Ccmtour  plot  of  the  Von  Mises  effective  stress  (psi)  in  HY-130 
steel  (a=1.537  in).  The  model  is  shown  as  2: 1  in  inches. 
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Figure  3.8b. 


Contour  i^ot  of  the  Von  Mises  effective  stress  (pi)  in  HY-130 
steel  (a=  1.623  in).  The  model  is  shown  as  2: 1  in  inches. 
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Figure  3.8c. 
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Contour  plot  of  the  Von  Mises  effective  stress  (psi)  in  HY-130 
steel  (a=1.7685  in).  The  model  is  shown  as  2: 1  in  inches. 
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Variation  of  7 along  the  crack  tip  ligament  for  7075  aluminum 
(0=1.537  in).  Note  the  variation  of  J  is  small  until  the  load  is  increased 
beyond  P*.  Note  that  1  in*lb/in''2  =  0. 17513  kJ/m''2. 
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Figure.  3. 10.  Variation  of  J  along  the  crack  tip  ligament  for  7075  aluminum 

(0=1.623  in).  Note  the  variation  of  7  is  small  until  the  load  is  increased 
beyond  P*.  Note  that  1  in»lb/in'^2  =  0.17513  kJ/m''2. 


0.312  P* 


Ligament  ratio  x/a  (in/in) 


Rgure  3.11.  Variation  of  J along  the  crack  tip  ligament  for  HY - 130  steel 
(<7=1.623  in).  Note  the  variation  of  J  is  small  until  the  load  i 
beyond  F*.  Note  that  1  in»lb/in^2  =  0.17513  kJ/m^2. 
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Figure  3. 12.  Variation  of  7  along  the  crack  tip  ligament  for  HY-130  steel 

(tf=1.7685  in).  Note  the  variation  of  7  is  small  until  the  load  is  increased 
beyond  P*.  Note  that  1  in*lb/in'^2  =  0. 17513  kJ/m'^2. 
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Metal  Matrix  Composites 

SCS-6/Ti-15-3,  was  chosen  as  our  model  system  due  to  the  technological 
significance  and  substantial  amount  of  published  experimental  work  on  this  material 
[48].  The  goal  in  this  segment  of  our  investigation  was  to  document  the  mechanics  of 
failure  under  various  interface  strengths.  As  discussed  earlier  in  Chapter  1,  certain 
failure  criteria  will  be  documented  and  utilized  to  determine  the  direction  and  occurrence 
of  material  failure.  The  local-global  analysis  of  SCS-6/Ti  15-3  w'as  studied  using  the 
ANSYS  finite  element  package. 

Macro-Mechanics  Analysis 


In  the  analysis  of  orthotropic  materials,  an  optimized  mesh  design  is  essential, 
even  more  so  than  in  the  well  documented  cases  of  steel  and  aluminum.  Four  mesh 
designs  with  126  to  456  elements  were  loaded  with  identical  boundary  conditions  as 
was  described  in  Chapter  2.  In  each  case,  the  effective  principal  stress  at  a  specified 
node  (r  =  3.74  mm,  0  =  8. 1*)  was  plotted  as  a  function  of  the  number  of  elements 
within  the  mesh  (see  Fig.  3.14).  This  point  is  far  enough  away  from  the  crack  tip  to 
provide  a  stable  solution.  As  a  further  check,  the  contour  plots  of  the  principal  stresses 
were  recorded  for  the  ^vo  meshes  with  the  highest  number  of  elements  (refer  to  Figs. 
3.15  and  3. 16).  Note  that  the  255  element  model  gives  drastically  different  results  from 
that  of  the  456  element  model.  Only  the  results  from  the  456  element  mesh  design  were 
used  in  the  boundary  condition  calculations 

To  accurately  model  the  experimentally  obsen'ed  failure  mechanisms,  a 
combination  of  failure  criteria  must  be  monitored.  The  initial  mechanisms  include  fiber 
pull-out  and  matrix  plasticity.  When  loads  are  further  increased,  the  fiber  or  matrix  may 
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Figure  3.15.  Effective  stress  contour  plot  in  255  element  macroscopic  model.  The 
plate  is  modeled  using  quaner  symmetry,  the  crack  tip  is  located  in  the 
lower  left  comer. 
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Figure  3. 16.  Effective  stress  contour  plot  in  456  element  macroscopic  model.  The 
plate  is  modeled  using  quarter  symmetr>',  the  crack  tip  is  located  in  the 
lower  left  comer. 
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fail  according  to  equations  1.22  and  1.23  respectively  as  was  described  in  Chapter  1. 
Since  fiber  pull-out  and  matrix  plasticity  may  precede  failure,  the  normal  stress  ratio  and 
Von  Mises  yield  criteria  will  be  used  to  predict  fiber  debonding  and  matrix  yielding 
respectively.  The  normal  stress  ratio  in  this  case  is  equal  to  the  transverse  stress  in  the 
composite  divided  by  a/,  the  tensile  strength  of  the  interface.  Thus  when  the  normal 
stress  ratio  exceeds  1.0  then  interfacial  failure  will  occur.  Experimental  results  [48] 
indicate  that  fiber  debonding  occurs  at  transverse  stresses  as  low  as  138  MPa,  well 
below  the  level  required  for  matrix  yielding  and  fiber  failure.  Thus  the  global  model 
took  into  consideration  material  nonlinearities,  but  the  solution  did  not  reach  the  required 
stress  state  for  plastic  deformation. 

As  was  described  in  detail  in  Chapter  2,  the  boundary  conditicms  for  the  micro¬ 
model  were  extracted  from  the  sdution  to  the  global  model.  The  critical  external  load 
required  to  initiation  fiber-matrix  debonding  was  reached  when  cfr=  138  MPa.  Note 
from  Table  2.2  that  this  is  well  below  the  stress  level  necessary  to  induce  yielding.  The 
induced  displacements  parallel  (UX)  and  perpendicular  (UY)  to  the  direction  of  the  crack 
were  then  plotted  and  fit  with  the  best  curve  as  shown  in  Figs  3. 17  and  3. 18.  The  four 
curves  (A,  B,  C  and  D)  represent  the  nodes  of  the  macro-model  as  was  depicted  in  Fig. 
2.7.  For  the  displacement  along  the  crack  tip  direction,  the  fitted  curves  are  given  in 
equations  3. 1  and  3.2  for  curves  C  and  D  respectively. 


y  =  -0.506  -  l.(M6xlO-4r  -  1.804x10-^ 

R2  =  0.972 

(3.1) 

y  =  -0.449  -  1.342x10^  -  1.627xlO-V 

R2  =  0.995 

(3.2) 

For  curves  A  and  B,  the  correlation  was  less  than  0.92,  therefore  \'isual  interpolation 
was  used  to  obtain  the  necessary  points.  For  the  displacements  normal  to  the  crack  tip 
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direction,  the  fitted  curv'es  are  given  in  equations  3.3-3.6  for  curves  A-D  respectively. 


y=  1.091x10-3 -5. 144xl0-3r- 

2.435xl0-5r2  -  4.987xl0-8r3  R2  =  0.999  (3.3) 


y  =  0.172  -  1.252xl0-3r  -  3.690xl0-6r2  - 


6.318xl0-9r3 . 3.506xl0-Hr4 

R2  =  0.999 

(3.4) 

y  =  0.253  -  8.688xl0-4r  -  1.269xl0-6r2 

R2  =  0.998 

(3.5) 

y  =  0.318  -  7.389x10^  -  7.572xl0-7r2 

R2  =  0.999 

(3.6) 

Micro-Mechanics  Analysis 

The  fiber,  matrix  and  interphase  regions  are  modeled  as  distinct  materials  on  the 

local  level.  Boundary  conditions  were  extracted  from  elements  of  the  global  model,  as 

was  described  in  the  previous  section,  whose  transverse  stresses,  o^,  exceed  the  stress 

required  to  initiate  fiber  debonding.  The  height  and  sides  of  the  local  model  are  each 
380  f*m.  Three  different  values  of  were  used  under  identical  boundary  conditions. 
Figure  3. 18  shows  the  deformed  mesh  when  oj  =  138  MPa.  The  outlined  area  indicates 
the  plane  of  the  crack.  Note  that  complete  fiber  debonding  occurred  along  all  three 
fibers.  When  the  strength  of  the  interface  was  increased  to  that  of  the  fiber,  extremely 
high  stresses  were  developed  within  the  fibers.  This  case  does  not  merit  additional 
discussion.  As  an  intermediate  value,  a/  was  set  equal  to  the  strength  of  the  matrix,  and 
load  is  effectively  transferred  to  the  fibers.  Figure  3.20  shows  the  stress  contours  when 
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Bgure  3. 17.  Displacements  obtained  from  the  solution  of  the  global  model.  The  data 
are  interpolated  using  polynomial  fits  on  layers  D  and  C,  while  visual 
interpolation  is  used  on  layers  B  and  A  (refer  to  Fig.  2.7). 


Displacement  parallel  to  the  crack  direction 
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Figure  3. 18.  Displacements  obtained  from  the  solution  of  the  global  model.  The  data 
are  interpolated  using  polynomial  fits  on  all  layers  (refer  to  Fig.  2.7). 
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Figure  3. 19.  Deformed  finite  element  mesh  showing  fiber  debonding  in 
SCS-6/Ti-15-3.  Crack  plane  is  outlin^  as  shown. 
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Figure  3.20.  Effective  stress  contour  in  SCS-6/Ti-15-3  when  Oi  =  934  MPa.  Note 
that  the  stresses  are  concentrated  within  the  fiber  due  to  effective  load 
transfer. 
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a,  =  Om.  In  this  case  the  maximum  stress  within  the  fiber  far  exceeds  its  ultimate 
strength,  indicating  instantaneous  fracture  of  the  fibers  under  the  assumed  nature  of  the 
interface. 

When  relatively  low  loads  (still  well  within  the  elastic  limit)  are  developed  about 
the  crack  tip  in  the  global  composite,  substantial  damage  occurs  on  a  microscopic  scale. 
As  was  observed  in  monolithic  steel  and  aluminum,  the  U-end  of  stress  decaying  away 
from  the  crack  tip  is  maintained  as  is  shown  in  Fig.  3.21.  When  the  signature  of  the 
stress  field  in  the  matrix  near  the  crack  tip  is  compared  to  the  signature  in  monolithic 
metals,  similarities  are  found.  Note  that  at  a  finite  distance  from  the  crack  tip,  the  stress 
increases  then  decreases  once  more. 

By  taking  a  closer  look  at  the  effective  stresses  near  the  interface  region,  it  is 
possible  to  document  the  mechanics  of  load  transfer.  Rgure  3.22  shows  the  orientation 
of  a  sectioned  surface  of  the  micro-model,  including  the  three  fibers  labeled  A,  B  and  C. 
Figures  3.23, 3.24  and  3.25  exhibit  how  load  is  transferred  from  the  matrix  to  the  fibers 
using  the  right  hand  definite  orientation  described  above.  These  three  figures  represent 
different  degrees  of  interfacial  bonding:  138  MPa,  534  MPa  and  934  MPa.  The  load 
transfer  ratio  is  defined  in  equation  3.7. 

Load  transfer  ratio  =  (3.7) 

Matrix 

From  these  plots,  it  is  evident  that  on  the  average,  the  broken  fiber  does  not  transfer  load 
as  effectively  as  does  the  undamaged  fibers.  Also  note  that  load  transfer  is  not  uniform 
and  that  it  depends  upon  the  angular  orientation  and  distance  from  the  crack  front.  In 
general,  the  interface  regions  on  the  plane  of  the  crack  front  transfer  more  load  than  does 
the  plane  furthest  from  the  crack  front 


Load  Transfer  Ratio  Load  Transfer  Ratio 

[Otlber/Omalrixl  lonber/f’matrixl 


Figure  3.23. 
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Part  A 


A  Fiber  A 
□  Fiber  B 
O  Fiber  C 


PartB 


A  Fiber  A 
□  HberB 
O  Fiber  C 


Part  A  shows  effective  load  transfer  in  fibers  A,  B  and  C  as  a 
function  of  angular  orientation  on  the  crack  plane.  Part  B  shows 
effective  load  transfer  on  the  plane  furthest  from  the  crack  plane.  The 
interface  strength  is  934  MPa. 


Load  Transfer  Ratio  Transfer  Ratio 
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Part  A 


A  Fiber  A 
□  Fiber  B 
O  Fiber  C 


Parts 


A  Fiber  A 
□  Fiber  B 
O  Fiber  C 


Figure  3.24.  Part  A  shows  effective  load  transfer  in  fibers  A,  B  and  C  as  a 

function  of  angular  orientation  on  the  crack  plane.  Part  B  shows 
effective  load  transfer  on  the  plane  furthest  from  the  crack  plane.  The 
interface  strength  is  534  MPa. 
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Effective 
Stress  in 
the  Matrix 
(MPa) 


Interface  Strength  (MPa) 


Figure  3.26.  Optimization  of  the  interface  strength  considering  the  stress  in  the  matrix 
near  the  fiber-matrix  interface  (£,  =  92.39  GPa).  The  optimum  \'alue  of 

the  interface  strength,  {oi)opt ,  may  be  obtained  by  reading  the  interface 
strength  at  the  intercept  of  Uie  two  curv'es. 
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An  estimate  of  the  optimum  interface  strength,  ( Oi)opt,  can  be  made  by  picking  a 
single  point  in  the  matrix  near  the  interface.  The  point  chosen  corresponds  to  the  node  at 
which  the  greatest  effective  stress  is  observed.  Figure  3.26  shows  that  as  the  degree  of 
interfacial  bonding  is  increased  beyond  730  MPa,  the  stress  in  the  matrix  exceeds  739 
MPa,  the  stress  in  the  matrix  at  the  fracture  strain  of  the  SCS-6  fibers  (0.8%). 

Therefore,  when  fabricating  SCS-6/Ti-15-3  composites,  it  would  be  ad\’antageous  to 
increase  the  degree  of  interfacial  bonding  to  730  MPa.  Keep  in  mind  that  this  analysis 
uses  an  assumed  value  of  the  interface  stiffness  equal  to  that  of  the  matrix  stiffness. 


CHAPTER  4 


SUMMARY  AND  CONCLUSIONS 

A  series  of  materials  have  been  evaluated  numerically  for  their  fracture  behavior 
and  response  to  applied  loads.  Monolithic  HY- 130  steel,  HY-140  steel,  5083-0 
aluminum  and  7075  aluminum  have  been  studied  using  ARLPAPST,  a  nonlinear  finite 
element  package.  In  the  case  of  HY-130  steel,  7075  aluminum  and  5083-0  aluminum, 
the  strain  energy  distribution  allows  estimates  of  the  crack  initiation  loads  and  the  failure 
load  to  be  made.  The  effect  of  crack  tip  plasticity  and  crack  tip  unloading  on  the  J- 
Integral  have  been  established.  Values  of  the  /-Integral  were  examined  for  a  range  of 
crack  to  width  ratios  (0.533  ^  o/Ws  0.884)  leading  from  small  scale  to  extensive 
plasticity.  When  the  loads  were  increased  above  a  critical  load  (which  is  attributed  to 
crack  initiation),  values  of  the  /-Integral  deviated  from  the  averaged  value  across  the 
ligament  of  the  specimens,  up  to  13%. 

The  damage  mechanics  in  the  composite  were  studied  using  a  local-global 
analysis.  The  boundary  conditions  for  the  three-dimensional  micro-model  were 
extracted  from  the  two-dimensional  macro-model  through  an  innovative  concept  of 
coordinate  rotation  and  interpdation.  The  results  of  the  micro-mechanical  modeling  are 
valid  when  considering  only  the  non-dimensionalized  quantities  and  not  the  absolute 
numbers.  The  convergence  test  for  the  micro-model  still  needs  to  be  performed.  Within 
this  context,  extensive  deformation  and  interfacial  debonding  were  found  near  the  crack 
tip  of  the  c(»nposite  even  under  the  apjMication  of  relatively  low  loads.  These  numerical 
predictions  confirmed  the  published  experimental  observations  that  SCS-6/Ti-15-3 
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delaminates  readily  under  loading.  Extensive  fiber  debonding  was  observ  ed  when  the 
interfacial  strength  was  set  to  138  MPa.  Loads  were  more  effectively  transferred  when 
the  interface  strength  was  increased  above  138  MPa.  Thus  the  optimum  strength  of  the 
interface  under  the  assumed  conditions  of  the  model  is  between  534  and  934  MPa.  The 
signature  of  the  effective  stresses  in  the  micro-mechanical  model  of  the  composite 
signifies  the  crack  propagation  direction  as  in  the  case  of  matrix  materials. 


CHAPTER  5 


RECOMMENDATIONS  FOR  FUTURE  RESEARCH 

The  recommendations  for  future  work  in  the  area  of  local-global  finite  element 
analysis  will  be  itemized  below: 

•  An  iterative  load  stepping  solution  procedure  is  needed  to  determine  the  exact  point  of 
fiber  failure  and/or  fiber  matrix  debonding. 

•  More  values  of  the  interface  tensile  strength,  shear  strength  and  stiffness  need  to  be 
analyzed  to  optimize  the  point  at  which  delamination  is  favored  over  fiber  fracture. 

•  The  micromechanical  model  is  designed  to  easily  accept  different  crack  lengths.  The 
effect  of  crack  length  on  delamination  and  fiber  fracture  should  be  studied. 

•  The  case  of  crack  bridging  should  be  studied  by  simply  varying  the  boundary 


conditions  of  the  micro-model. 
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APPENDIX  A 


J-INTEGRAL  CUTTING  PROGRAM 


C  J-Integral  cutting  program  (JCUT.FOR)  created  by  Mike  House  on 
C  20-SEP-1990 11:19:07.02 

C 

To  cut  the  J-Integral  values  directly  from  the  main 
PAPST  output.  Program  can  handle  any  number  of 
increments. 


C  »>OBJECriVE«< 
C 
C 
C 

C  »>COMMENTS«< 
C 


Must  change  the  input  and  output  files  to  be  accessed 
by  program  within  the  open  statements. 


CHARACTER*!  A(132) 

CHARACTER  B»10.LOADS*40 
CHARACTER*  1 1JINT(  15,50), CHECK 

OPEN(UNIT=18,HLE='[mbh.monolith.HY-130.PAPST6}papst6.out', 

.STATUS='OLD') 

OPEN(UNIT=21,HLE='[mbh.monolith.HY-130.PAPST61jl30-3.OUT', 

.STATUS='NEW') 

REWIND  18 
WRITE(*,10) 

10  FORMAT(5X,'INPUr  NUMBER  OF  PATHS  CALCULATED  FOR  J- 

.INTEGRAL') 

READ  ».  NJINT 
WRITE(*,40) 

40  FORMAT(5X,lNPUT  NUMBER  OF  LOAD  INCREMENTS') 

READ  *,  NINC 
C - 

C  CUT  DATA  FROM  MAIN  OUTPUT 
C - 


DO  100  I=1,NINC 

DO  200  J=l, NJINT 

3  READ(  18,5000)  B,CHECK 

5000  FORMAT(T12,A10,T39,A11) 
IF(B.EQ.'J  INTEGRAL')  THEN 

JINT(J,I)=CHECK 
GO  TO  200 

ELSEIF(B.EQ. 'INCREMENT  ')  THEN 
DO  5003  JJ=1,NINC+1 

READ(  18,5001)  LOADS 

5001  FORMAT(T14,A40) 

WRITE(21,5002)  LOADS 
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5002  FORMAT(T10.A40) 

5003  CONTINUE 
GO  TO  3 

ELSE 

GO  TO  3 

ENDIF 

200  CONTINUE 

100  CONTINUE 

C - 

DO450JCOD=l.NINC 

WRITE(21,110)  JCOL 

1 10  FORMAT(/T7,’VALUES  OF  THE  J-INTEGRAL  FOR  I NCREMENT 

.'I3//T3.'PATH'.T10.'J-1NTEGRAL'/T3.' — '.TIO,' . 7) 

DO440IROW=l.NJINT 

WRITE(21.460)  IROW4INT(lROW.JCOL) 

460  FORMAT(T4,I2.T  10^  1 1) 

440  CONTINUE 

450  CONTINUE 
STOP 
END 


APPENDIX  B 


POSTSCRIPT  CONNECTIVITY  PROGRAM 


C  POSTSCRIPT  Connectivity  program 

C  created  by  Mike  Houston  20-AUG- 1991  19:52:11.02 

C 

C  »>OB  JECTIVE<<< 

C 

C  »>VARIABLES«< 

C 
C 
C 

C  »>COMMENTS«< 

C 
C 
C 

REAL*8X(50(X)).Y(5000).MAXSPW,DEFX(5000.50),DEFY(5000,50) 
CHARACTER*  10  JUNK^UNKS 
CHARACTER*  1 1  JUNK2 

INTEGER  C1(500),C2(500),C3(500),C4(500),MAT(500).DEFORM 
C 

READ(7,12)  NN,NE,NM 
12  FORMAT(/3I5) 

ISKIP  =  NM  *2 
DO  14 1=1, ISKIP 
READ(7,15)JUNK 
15  FORMAT(AIO) 

14  CONTINUE 
C 

MAXSPW  =  0.000001 
DO  1000 1=1, NN 

READ(11,20)X(I),Y(I) 

IF(X(I).GE.MAXSPW)  MAXSPW=X(I) 

20  FORMAT(T  17,F9.7, 1 1X,F9.7) 

1000  CONTINUE 

C  Code  now  calculates  size  of  grid  to  place  onto  paper 
C 

IWIDTH  =  432.0/MAXSPW 
C 

C  Postscript  recognizes  72.0  dpi,  therefore  there  are  72*8.5  columns 
C  432.0  allows  for  1.5"  left  margin  and  1.0"  right  margin 
C 


To  generate  model  meshes  showing  connectivity 

NN  =  Number  of  nodes  in  the  mesh 
NE  =  Number  of  elements  in  the  model 
MAXSPW  =  Specimen  width  in  inches 

Necessary  files  includes  the  raw  data  file,  all  comer 
nodes  and  their  connectivity  (this  is  usually  the  data 
submitted  to  PAPST),  a  file  containing  the 
coordinates  of  the  comer  nodes,  and  the  output  file. 
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WRITE(10.10)  IWIDTH.I WIDTH 
C 

WRITE(6,16) 

16  FORMAT(T5,'Is  data  in  2/comer  mode  or  full  mode?  (l=2/comer)') 
READ(5,»)  MODE 
IF(MODE.EQ.l)THEN 
D06000J=1,2 

READ(7,18)  C  1(J),C2(J).C3(J).C4(J),MAT(J) 

WRITE(10,2)  X(C1(J)).Y(C1(J)),X(C2(J)).Y(C2(J)),X{C3(J)), 
.Y(C3(J)).X(C4(J)).Y(C4(J)).X(C1(J)).Y(C1(J)) 

6000  CONTINUE 
D0  2000J=3,NE 

READ(7. 1 020)  C 1  (J),C2(J),C3(J),C4(J).MAT(J) 

WRITE(10,2)  X(C1(J)).Y(C1(J)).X(C2(J)).Y(C2(J)),X(C3(J)), 
.Y(C3(J)).X(C4(J)),Y(C4(J)).X(C1(J)),Y(C1(J)) 

2000  CONTINUE 

ELSEIF(MODE.NE.l)  THEN 
DO  19J=1.NE 

READ(7.18)  C  1(J).C2(J).C3(J),C4(J),MAT(J) 

WRITE(10.2)  X(C1(J)).Y(C1(J)),X(C2(J)),Y(C2(J)).X(C3(J)). 
.Y(C3(J)),X(C4(J)).Y(C4(J)).X(C1(J)).Y(C1(J)) 

19  CONTINUE 
ENDIF 
WRITEdO^ 

5  FORMAT(T2,'  0.6  setgray') 

C 

D046J=1.NE 

IF(MAT(J).EQ.2)  THEN 

WRITE(10,4)X(C1(J)).Y(C1(J)),X(C2(J)),Y(C2(J))^(C3(J)). 

.Y(C3(J)),X(C4(J)),Y(C4(J))^(C1(J)),Y(C1(J)) 

ENDIF 

46  CONTINUE 
C 

C  Now  we  want  to  read  in  displaced  nodal  values  and  add  to  X  and  Y 
C 

WRITE(6,4002) 

4002  FORMAT(T5,'FYint  deformed  mesh?  ( Isyes)*) 

READ(5  »)  DEFORM 
IF(DEFORM.EQ.l)  THEN 

GOTO  9998 

ELSE 

GOTO  9999 

ENDIF 

9998  WRITE(6,4003) 

4003  FORMAT(T5,'How  many  ccmverged  solutions?’) 

READ(5,*)  INUM 

DO  6001 11=1,NN 

DEFX(I1,1)=X(I1) 

DEFY(I1,1)=Y(I1) 

6001  CONTINUE 
C 
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D043(X)K=1.INUM 

READ(  15,4000)  JUNK2 
FORMAT(TlO,An) 

IF(JUNK2.EQ.'NODE  NUMBER')  THEN 
REACH  15.*) 

D0  4100J=1.NN 

REACH  15.4010)  DX,DY 
FORMATa25,El  1.8,T45,E1 1.8) 
DEFX(J,K)=X(J)+DX 
DEFY(J,K)=Y(J)+DY 
CONTINUE 
GOTO  4300 
ENDIF 
GOTO  4030 

4300  CONTINUE 
C 

C  Now  develop  format  for  datafile 
C 

4200  WRITE(6,4201) 

4201  FORMAT(T5,'Overlay  which  converged  solution?*) 

REACH5.*)  IC 

DO  4301  J=1,NE 

WRITE(  10,2)  DEFX(C1(J).1C),DEFY(C1(J),IC),DEFX(C2(J),IC), 
.DEFY(C2(J).IC).DEFX(a(J),IC),DEFY(C3(J),IC),DEFX(C4(J),IC), 
.DEFY(C4(J),IC).DEFX(C  1  (J).IC),DEFY(C  1(J),IC) 

4301  CONTINUE 

2  FORMAT(2X,F10.8,1X,F10.8.'  moveto  ',F10.8,1X,F10.8.’  lineto  '/ 

.1X,F10.8,1X,F10.8,'  lineto  ’,F10.8,1X.F10.8,’  lineto  V1X,F10.8. 

.  1X,F10.8,'  lineto  stroke') 

4  FORMAT(2X,F10.8,1X,F10.8,’  moveto  ’,F10.8,1X.F10.8,’  lineto  '/ 

.1X,F10.8,1X,F10.8,'  lineto  ',F10.8,1X,F10.8,'  lineto  '/1X,F10.8. 
.1X,F10.8,'  lineto  fill') 

C 

C  This  is  for  all  nodes  given 
C 

18  FORMAT(T8,I3.T23,I3,T38,I3,T53,I3,T70,1 1 ) 

10  FORMATC  %!'/'  103  320  translate’/’  0  seUinewidth’/ 

.T2.I4,lx,I4'  scale') 

C 

C  This  is  for  only  comer  nodes  given 
C 

1020  FORMAT(T8,I3,T13.I3.T18,I3,T23,I3.T70,1 1) 

9999  WRITE(103) 

3  FORMATC  showpage') 

STOP 

END 


4030 

4000 

4010 

4100 


APPENDIX  C 


PROGRAM  LISTING  OF  2-D  ORTHOTROPIC  PROBLEM 


/SHOW,Xll 

/PREP7 

/title,  2-D  MMC  Laminate 

CSYS.O 

KAN.O 

KAY. 6,0 

ET.1.82...3.,2 

R,THICKNESS.1.55E-3 

K.  1,0,0 

K.2.2.65e-3.0 

K.3.3.175e-3.0 

K,4,3.7e-3,0 

K.5.6.35e.3,0 

K.6.0.0.525e-3 

K,7,2.65e-3,0.525e.3 

K,83.175e-3,0.525e-3 

K.9,3.7e-3,0.525e-3 

K,10,6.35e-3,0.525e-3 

K,ll,0.1.525e-3 

K,12,2.65e-3,1.525e-3 

K,13,3.175e.3,1.525e-3 

K,143.7e-3,1.525e-3 

K.15,6.35e-3,1.525e-3 

K.16,0,12.7e-3 

K.I7,2.65e-3,12.7e.3 

K,18,3.175e-3,12.7e-3 

K,19,3.7e-3,12.7e-3 

K, 20,6.35e-3,12.7e-3 
/VIEW.l. 1,1,1 

L. 2.1.3 
L,2.3,3 
L.3,4,4 
L.4.5.3 
L,7,6,3 
L,7.8.3 
L.8,9.4 
L.9.103 
L.12.11,3 
L.12,13,3 
L,13.14,4 


!  Define  cartesian  coordinates 
!  Define  static  analysis 
!  Small  deflection  option 


!  Define  keypdnts  1-8 


!  Change  perspective 
!  List  nodes  as  from,to,div,ratio 
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L,  14. 15,3 
L, 17, 16,3 
L,  17, 18,3 
L.  18, 19,4 
L.  19,20.3 

L,  1 ,6,3  !  Up  and  down  segments 

L,2,7,3 

L,3,8.3 

L.4,9,3 

L,5,10,3 

L.6,11.2 

L,7,12,2 

L.8,13.2 

L,9.14,2 

L,10.15.2 

L,1 1,16,10,3 

L,12,17.10,3 

L.13, 18,10,3 

L.  14, 19. 10,3 

L.  15,20, 10,3 

ELS1ZE„2  ’.Default  number  of  elem  divisions 

A,1.2,7,6 

A,2,3,8,7 

A,3,4,9,8 

A.4,5,10,9 

AMESH,1,4 

ELSIZE„2 

A,6,7,12,ll 

A,7,8,13.12 

A,8,9,14,13 

A,9,10,15,14 

AMESH,5,8  !  MESH  AREAS  5  TO  8 

A,ll,12,17,16 
A,12,13,18,17 
A,13,14,19,18 
A,  14, 15,20, 19 
AMESH,9,12 
ITER,1,1.1 
KBC  0 

MP.EX,1,122.7266E9 
MP,EY,1,208.222E9 
MP.NUXY,  1,0.299 
LSBC.22 
LSBC,17 
LSBC.27 
LSBC,3 
LSBC,4 

KD,  16.UY.2.4983954e-6„  1 
KD.  1 7,UY  ,2.4983954e-6,.  1 
KD,  18,UY  ,2.4983954e-6„  1 
KD,  19.UY,2.4983954e-6„  1 


!  DEHNES  SYMMETRIC  BC'S 
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KD,20,U  Y  .2.4983954e-6..  1 

SBCTRA  !  TRANSFERS  BC'S  TO  MESH 

WAVE 


